The non-simplified border is 3229 points, while the simplified border is 487 points. One consequence of doing this computationally is you do not get to choose where or what is simplified.
| Tesselation | # of Features | Mean Area (km^2) | Std. Deviation | Total Area (km^2) |
|---|---|---|---|---|
| Normal County | 3108 | 2521.745 | 3404.325 | 7837583 |
| Voronoi Tesselation | 3107 | 2521.830 | 2886.800 | 7835326 |
| Triangulation Tesselation | 6196 | 1251.865 | 1576.180 | 7756553 |
| Grid | 3108 | 2728.126 | 0.000 | 8479014 |
| Hex Grid | 2271 | 3763.052 | 0.000 | 8545891 |
Different tessellations seemed to either increase or decrease the density of dams. For example, the Voronoi plot seemed to put more density on some northwestern counties, where the triangulated and hex grid do not have as much density in the same locations as comparison. This is related to the MAUP because as the size and coverage of each tessellation changes, they may generalize parts of the data, resulting in bias. I will choose the hex grid because it has the smallest # of features and generalizes the data the least.
I am choosing Flood Control, Fire Protection, Water Supply, and Irrigation dams because I feel these are some of the highest prority dams in terms of their uses.
As we can see, the flood control dams seem to be focused along the Mississippi river, which comes to no surprise. The fire protection da,s also seem to be in the same location, but more focused. The water supply dams are in the north-central location of the US. and the irrigation dams are located in two hot spots, central USA and the East Coast.
knitr::opts_chunk$set(echo = FALSE, warning = F, message = F)
options(dplyr.summarise.inform=F)
library(sf)
library(tidyverse)
library(leaflet)
library(leafpop)
library(USAboundaries)
library(rmapshaper)
library(mapview)
library(gridExtra)
library(units)
library(knitr)
library(readxl)
library(gghighlight)
source('../R/utils.R')
counties <- us_counties() %>% st_transform(5070) %>% filter(!(state_name %in% c('Alaska', 'Puerto Rico', 'Hawaii')))
cent_county <- st_centroid(counties)
cent_county <- cent_county %>% st_combine()
# Voronoi
voronoi <- cent_county %>% st_voronoi() %>% st_cast() %>% st_as_sf() %>% mutate(id = 1:n())
# Triangulate
triangulate <- cent_county %>% st_triangulate() %>% st_cast() %>% st_as_sf() %>% mutate(id = 1:n())
# Square Grid
grid <- counties %>% st_make_grid(n = 70) %>% st_cast() %>% st_as_sf() %>% mutate(id = 1:n())
# Hexagonal Grid
hex_grid <- counties %>% st_make_grid(square = FALSE, n = 70) %>% st_cast() %>% st_as_sf() %>% mutate(id = 1:n())
border_raw <- counties %>% st_union()
border <- ms_simplify(border_raw, keep = .15)
plot(border)
#npts(border_raw)
#npts(border)
# Voronoi
voronoi <- st_intersection(voronoi, border)
# Triangulate
triangulate <- st_intersection(triangulate, border)
plot_map <- function(sf_object, title ){
ggplot() +
geom_sf(data = sf_object, fill = 'white', col = 'navy', size = .2) +
theme_void() +
labs(title = paste0(title, ' Map'),
caption = paste0('There are ', count(sf_object), ' features in the tesselation'))
}
counties_map <- plot_map(counties, 'Normal County')
voronoi_plot <- plot_map(voronoi, 'Voronoi Tesselation')
triang_plot <- plot_map(triangulate, 'Triangulation Tesselation')
grid_plot <- plot_map(grid, 'Grid')
hexgrid_plot <- plot_map(hex_grid, 'Hex Grid')
grid.arrange(counties_map, voronoi_plot, triang_plot, grid_plot, hexgrid_plot, nrow = 3)
tess_summary <- function(sf_object, tess){
area <- st_area(sf_object) %>% set_units('km^2') %>% drop_units()
num_feat <- length(area)
total_area <- sum(area)
avg_area <- mean(area)
std_area <- sd(area)
summary<- data.frame(Tesselation=tess,
`# of Features` = num_feat,
`Mean Area (km^2)` = avg_area,
`Std. Deviation` = std_area,
`Total Area` = total_area)
return(summary)
}
summary_table <- rbind(
tess_summary(counties, 'Normal County'),
tess_summary(voronoi, 'Voronoi Tesselation'),
tess_summary(triangulate, 'Triangulation Tesselation'),
tess_summary(grid, 'Grid'),
tess_summary(hex_grid, 'Hex Grid')
)
kable(summary_table, col.names = c('Tesselation', '# of Features', 'Mean Area (km^2)', 'Std. Deviation', 'Total Area (km^2)' ),
caption = 'Tesselation Summary Table')
conus <- state.abb[!(state.abb %in% c('AK', 'HI', 'PR'))]
dams_raw <- read_excel('../data/NID2019_U.xlsx') %>% filter(STATE %in% conus) %>% filter(LONGITUDE != 0 | LATITUDE != 0)
dams <- dams_raw %>% filter(!is.na(LONGITUDE) & !is.na(LATITUDE)) %>% st_as_sf(coords = c('LONGITUDE','LATITUDE'), crs = 4326) %>% st_transform(5070)
point_in_polygon = function(points, polygon, id){
st_join(points, polygon) %>%
st_drop_geometry() %>%
count(.data[[id]]) %>%
setNames(c(id, "n")) %>%
left_join(polygon, by = id) %>%
st_as_sf()
}
plot_pip = function(data, title){
ggplot() +
geom_sf(data = data, aes(fill = n), alpha = .9, size = .2, col = NA) +
scale_fill_viridis_c() +
theme_void() +
theme(plot.title = element_text(face = "bold", color = "navy", hjust = .5, size = 24)) +
labs(title = title,
caption = paste0("Number of Dams: ", sum(data$n) ))
}
dams_normal <- point_in_polygon(dams, counties, 'geoid')
dams_voronoi <- point_in_polygon(dams, voronoi, 'id')
dams_triang <- point_in_polygon(dams, triangulate, 'id')
dams_grid <- point_in_polygon(dams, grid, 'id')
dams_hex <- point_in_polygon(dams, hex_grid, 'id')
county_dam_map <- plot_pip(dams_normal, 'Dams per County')
voronoi_dam_map <- plot_pip(dams_voronoi, 'Dams per Voronoi Tile')
triang_dam_map <- plot_pip(dams_triang, 'Dams per Triangle Tile')
grid_dam_map <- plot_pip(dams_grid, 'Dams per Grid')
hex_dam_map <- plot_pip(dams_hex, 'Dams per Hex-Grid')
grid.arrange(county_dam_map, voronoi_dam_map, triang_dam_map, grid_dam_map, hex_dam_map, nrow = 3)
flood_dams <- dams %>% filter(grepl('C', PURPOSES))
fire_dams <- dams %>% filter(grepl('P', PURPOSES))
supply_dams <- dams %>% filter(grepl('S', PURPOSES))
irri_dams <- dams %>% filter(grepl('I', PURPOSES))
pip_flood <- point_in_polygon(flood_dams, hex_grid, 'id')
pip_fire <- point_in_polygon(fire_dams, hex_grid, 'id')
pip_supply <- point_in_polygon(supply_dams, hex_grid, 'id')
pip_irri <- point_in_polygon(irri_dams, hex_grid, 'id')
flood_map <- plot_pip(pip_flood, 'Flood Control Dams') + gghighlight(n > (mean(n) + sd(n)))
fire_map <- plot_pip(pip_fire, 'Fire Protection Dams') + gghighlight(n > (mean(n) + sd(n)))
supply_map <- plot_pip(pip_supply, 'Water Supply Dams') + gghighlight(n > (mean(n) + sd(n)))
irri_map <- plot_pip(pip_irri, 'Irrigation Dams') + gghighlight(n > (mean(n) + sd(n)))
grid.arrange(flood_map, fire_map, supply_map, irri_map, nrow= 2)
rivers <- read_sf('../data/majorrivers_0_0/MajorRivers.shp')
miss <- rivers %>% filter(SYSTEM == 'Mississippi')
top_dams <- dams %>% filter(HAZARD == 'H') %>% group_by(STATE) %>% slice(which.max(NID_STORAGE)) %>% ungroup() %>%
st_transform(4326)
leaflet() %>%
addProviderTiles(providers$OpenStreetMap) %>%
addCircleMarkers(data = top_dams, color = NA,
fillColor = 'red',
radius = top_dams$NID_STORAGE/1500000,
popup = popupTable({top_dams %>% st_drop_geometry() %>%
select(State = STATE,
`Dam Name` = DAM_NAME,
Storage = NID_STORAGE,
Purposes = PURPOSES,
`Year Completed` = YEAR_COMPLETED)},
feature.id = F, row.numbers = F)) %>%
addPolylines(data = miss)